function [sed_output] = stat_ensemble_dynamics(data,win,bin,null,splitmin,stepsize,alpha,berr,sorted)
% This code calculates both the Bayes's estimator for KL-divergence, according
% to a user specified null hypothesis, and the population firing rate of the
% ensemble. For a more detailed description of the Bayes's estimator see
% 'kdq_bayes_KL.m'. Also, check 'population_fr.m' for more details on the
% calculation of the population firing rate.
%
% This is code from the paper "A Statistical Description of Neural 
% Ensemble Dynamics," Long and Carmena 2011, submitted to  
% Frontiers in Computational Neuroscience.
% by: John D. Long II
% contact: jlong29@berkeley.edu or jlong29@gmail.com.
% Download at: http://code.google.com/p/kdq-bayes-kl
%
% INPUT:
% data - an MxN matrix of M time points and N units
% win  - an integer specifying the number of samples to include in the sliding
%        window.
% splitmin - an integer specifying the number of samples to include in a bin 
%       before splitting the bin via the kdq-tree (see 'Bern_tree.m').
% stepsize - an integer specifying the number of samples to slide between each
%       window evaluations.
% alpha - a real number setting the Dirichlet prior parameters to the same
%         value. Default: alpha = 1 (uniform prior). alpha < 1 biases output
%         toward non-uniform posterior and alpha > 1 biases output toward
%         uniform posterior.
% null - a string indicating which null hypothesis to evaluate (options below)
% berr - if this input is not empty, the standard deviation of the Bayes's
%        estimator will be calculated at each point. This is computationally 
%        expensive.
% 
% null hypotheses: 
% fixed - deviation of subsequent samples from an initial window's sample
%         distribution, syntax: 'fixed'
% derivative - change in KL-Divergence between adjacent windows slid by
%              stepsize, syntax: 'derivative','dkl' (default)
% independent - compares each sample within a window against a surrogate 
%               independent sample--generated by permuating each channel's 
%               time indices to break correlations--to test for correlations in
%               the data, syntax: 'independent','ind'
% different correlations - tests adjacent sliding windows for differences in the
%                          correlation structures. This test is 2 nested
%                          hypotheses:
%                          1) Are the samples the 2 adjacent windows 
%                             independent or correlated?
%                          2) Are the samples in adjacent windows from the same
%                             distribution?
%                          syntax, 'dcorr'
%
% OUTPUT
% The output of this code is formatted to be compatible with:
% 'stat_ensemble_dynamcis_visualization.m'
% sed_output is a stucture with the following fields
%   .data - a compact version of the original data, used for later visualization
%   .T - The number of data points
%   .N - The number of variables
%   .win - the win parameter
%   .bin - the bin parameter
%   .KLs - The KL-divergence estimated
%   .m_KL - the mode of the time-series of KL-divergence values
%   .std_KL - the standard deviation of the time-series of KL-divergence values
%   .pop_fr - the population firing rate estimated
%   .m_pop_fr - the mean of the time-series of population firing rate values
%   .std_pop_fr - the standard deviation of the time-series of pop_fr

% Argument check
[dp,N] = size(data);
if dp < N
    error(['More variables than datapoints. Check tranpose: data should be M '...
    'datapoints by N variables'])
end
if nargin < 2 || isempty(win)
    win      = 200;
    disp(['Win is set at ' num2str(win) '.'])
end
if nargin < 3 || isempty(bin)
   bin = 20;
    disp('Bin set at 20 samples.')
end
if nargin < 4 || isempty(null)
    null = 'dkl';
    disp('No null hypothesis selected: Evaluating KL-Divergence between sliding windows.')
end
if nargin < 5 || isempty(splitmin)
    splitmin = 5;
    disp(['Splitmin is set at ' num2str(splitmin) '.'])
end
if nargin < 6 || isempty(stepsize)
    stepsize = 1;
    disp(['Stepsize is set at ' num2str(stepsize) '.'])
end
if nargin < 7 || isempty(alpha)
    alpha = 0.5;
    disp(['Dirichlet alpha set at ' num2str(alpha) ': Jeffrey''s prior.'])
end
if nargin < 8 || isempty(berr)
    berr = [];
    disp('The standard deviation of the bayes''s estimator will not be calculated.')
end
if nargin < 9
    sorted = 'sort';
end
if (dp/bin) < win
    error(['Your window is bigger than the number of available data points. ' ...
        'There is no analysis to run.'])
end


% check whether to store sorted data or not
if ~isempty(sorted)
    fprintf('Storing data sorted by firing rate.\n\n')
    [temp,ind] = sort(sum(data));
    data = data(:,ind);
    % store efficient representation of spike ensemble data
    [rind,cind] = find(data);
else
    fprintf('Storing unsorted data.\n\n')
    % store efficient representation of spike ensemble data
    [rind,cind] = find(data);
    [temp,ind] = sort(sum(data));
    data = data(:,ind);
end
    
    
%%% KL-divergence %%%
%%%%%%%%%%%%%%%%%%%%%
if bin ~=1
    [data1]     = bin_spikes(data,bin,1);
else
    data1 = data;
end
[KLs, null] = kdq_bayes_KL(data1,win,splitmin,stepsize,alpha,null,berr);
% Store output stats
if strcmp(null.model,'different correlations')
    KLs1    = KLs.ind;
    KLs2    = KLs.dkl;
    md_KL1  = null.mode_ind;
    mn_KL1  = null.mean_ind;
    md_KL2  = null.mode_dkl;
    mn_KL2  = null.mean_dkl;
    std_KL1 = null.std_ind;
    std_KL2 = null.std_dkl;
else
    KLs1    = KLs;
    KLs2    = NaN;
    md_KL1  = null.mode;
    mn_KL1  = null.mean;
    md_KL2  = NaN;
    mn_KL2  = NaN;
    std_KL1 = null.std;
    std_KL2 = NaN;
end

% Population Firing Rate
% counts for population firing rate
if bin ~=1
    [data2] = bin_spikes(data,bin,0);
else
    data2 = data;
end

%%% Population firing rate %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pop_fr = population_fr(data2,win,bin);

% Estimate variance of this estimate
% Using shuffled data
ind       = randsample(1:size(data,1),size(data,1));
shuff     = data2(ind,:);
pop_fr_rs = population_fr(shuff,win,bin);

% Resample mean and standard error
m_pop_fr   = mean(pop_fr_rs(win:end));
std_pop_fr = std(pop_fr_rs(win:end));

% SET OUTPUT FIELDS
% General data information
sed_output.data       = [rind cind];
sed_output.T          = size(data,1);
sed_output.N          = size(data,2);
sed_output.win        = win;
sed_output.bin        = bin;
% KL-divergence information
sed_output.KLs1       = KLs1;
sed_output.KLs2       = KLs2;
sed_output.md_KL1     = md_KL1;
sed_output.mn_KL1     = mn_KL1;
sed_output.md_KL2     = md_KL2;
sed_output.mn_KL2     = mn_KL2;
sed_output.std_KL1    = std_KL1;
sed_output.std_KL2    = std_KL2;
sed_output.null_model = null.model;
% Population firing rate information
sed_output.pop_fr     = pop_fr;
sed_output.m_pop_fr   = m_pop_fr;
sed_output.std_pop_fr = std_pop_fr;